home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Power Programmierung
/
Power-Programmierung (Tewi)(1994).iso
/
magazine
/
drdobbs
/
1991
/
05
/
co_dsp.asc
< prev
next >
Wrap
Text File
|
1991-03-15
|
21KB
|
529 lines
_ADDING THE POWER OF DSP TO YOUR APPLICATIONS_
by Jim Bittman
[LISTING ONE]
/*-------------------- DDJ.C-----------------------------------------
| PC-Resident program for FFT and for controlling DSP co-processor
| The Makefile for this program is:
| .c.obj:
| cl /AL -c -DANSI $<
| ddj.exe: ddj.obj dspif.obj dspif.h four1.obj realft.obj
| cl /AL ddj.obj dspif.obj four1.obj realft.obj \
| -link dsputill.lib hqmcil.lib
| This program assumes the DspHq environment and requires additional
| support files (menu definition and function specification files)
| not listed in the Makefile. However, it should not be difficult to
| excerpt the relevant portions of this code for standalone execution.
+------------------------------------------------------------------*/
#include <stdio.h> /* include for NULL define */
#include <math.h> /* include for math functions */
#include "hq_ci_.h" /* include for dsphq interface */
#include "nr.h" /* include nr function prototypes */
#include "nrutil.h" /* include nr utility prototypes */
#include "dsputil.h" /* include CAC prototypes */
#include "dspif.h" /* dsp interface function header */
/*---------- define the menu parameter types ----------------------*/
typedef unsigned long parm1_type; /* Buffer 1 */
typedef unsigned long parm2_type; /* Buffer 3 */
typedef menu_float parm3_type; /* Cosine Amplitude */
typedef menu_float parm4_type; /* Cosine DC-offset */
typedef menu_float parm5_type; /* Cosine Frequency */
typedef menu_float parm6_type; /* Cosine Sammple Rate */
/*------------ constant definitions -------------------------------*/
#define SQR(a) ((a)*(a))
#define PI (float) 3.141592654
#define BAD_FUNC_NUM (int) -1 /* user error code */
/*----- Function Constants For Dsp32 'C' & Assembly Code ------*/
#define DSP32_X 1 /* converts IEEE ==> DSP numbers */
#define IEEE32_X 2 /* converts DSP ==> IEEE numbers */
/*----- Function Constants For Dsp32 'C' Code ----------------*/
#define GENCOS_C 3 /* generate test cosine */
#define REALFT_C 4 /* performs realft, from numerical recipes */
#define IREALFT_C 5 /* performs inv_realft, from numerical recipes */
#define LOGMAG_C 6 /* calculates log-magnitude */
#define RFFTA_C 7 /* calculates real fft using app-lib */
/*----- Function Constants For Dsp32 Assembly Code ------------*/
#define RFFT_S 3 /* performs rfft, from AT&T application library */
#define MAG_S 4 /* calculates magnitude-squared */
#define LOG_S 5 /* calculates log */
/*----- Function Prototypes ----------------------------------*/
int ilog2 (unsigned int);
void log_mag (float far * i1, float far * o1, long bs);
void scale_data (float far *output1, float scale_val, long np);
void synth_cos (float far *data1, long np, float a, float d, float f,
float s);
main(int argc, char *argv[])
{ long indx; /* used for loop index */
int func_num; /* for function number from dsphq */
long np; /* for block size from dsphq */
float far * input1; /* array address */
float far * output1; /* array address */
long bs; /* address of DSP blocksize var */
long flag; /* address of DSP funcnum flag */
parm1_type b1; /* DSP Buffer #1 */
parm2_type b2; /* DSP Buffer #2 */
parm3_type amp; /* cosine amplitude */
parm4_type dco; /* cosine DC offset */
parm5_type freq; /* cosine frequency */
parm6_type samprate; /* cosine sample rate */
init_intfc(argc, argv); /* init dsphq interface */
func_num = get_func_num(); /* get the function number */
np = get_block_size(); /* get the block size */
/* get menu parameters */
b1 = get_parm(1); /* DSP Buffer #1 */
b2 = get_parm(2); /* DSP Buffer #2 */
amp = get_parm(3); /* cosine amplitude */
dco = get_parm(4); /* cosine DC offset */
freq = get_parm(5); /* cosine frequency */
samprate = get_parm(6); /* cosine sample rate */
/* get array addresses */
input1 = get_data_in_ptr(1); /* base address of input #1 */
output1 = get_data_out_ptr(1); /* base address of output #1 */
/* perform selected function */
switch (func_num)
{
case 1 : /*--- Synthesize Cosine Using PC ---------------------*/
synth_cos(output1, np, amp, dco, freq, samprate);
break;
case 2 : /*---- Numerical Recipes Forward Real FFT Using PC ----*/
output1--; /* NR funcs index at 1 */
realft(output1, np>>1, 1); /* forward real fft */
break;
case 3 : /*---- Numerical Recipes Inverse Real FFT Using PC ----*/
output1--; /* NR funcs index at 1 */
realft(output1, np>>1, -1); /* inverse real fft */
output1++; /* reallign address */
scale_data(output1,1.0/(np >> 1),np); /* restore original ampl. */
break;
case 4 : /*------ Calculate LOG(10)-[MAGNITUDE] using PC --------*/
if (input1)
log_mag(input1, output1, np); /* take logmag of input */
else
log_mag(output1, output1, np); /* perform logmag in-place */
break;
case 5 : /*--------- Synthesize Cosine Using DSP-32C -------------*/
init_dsp("ddj_32c.32c",&flag,&bs,np); /* download dsp code & init */
dsp_dl_fp(get_addr("amp"),amp); /* download floats */
dsp_dl_fp(get_addr("dco"),dco);
dsp_dl_fp(get_addr("freq"),freq);
dsp_dl_fp(get_addr("samprate"),samprate);
set_dspbuf("o1", b1); /* set dsp buffer address */
dsp_dl_int(flag,GENCOS_C); /* invoke function on dsp */
wait4dsp(flag); /* wait for dsp to finish */
ul_float(output1,np,flag,b1); /* upload results */
break;
case 6 : /*---- Numerical Recipes Forward Real FFT Using DSP-32C ----*/
init_dsp("ddj_32c.32c",&flag,&bs,np); /* download dsp code & init */
set_dspbuf("o1", b1); /* set dsp buffer address */
dl_float(input1,np,flag,b1); /* download float array */
dsp_dl_int(flag,REALFT_C); /* execute "realft" on dsp */
wait4dsp(flag); /* wait for dsp to finish */
ul_float(output1,np,flag,b1); /* upload results */
break;
case 7 : /*---- Numerical Recipes Inverse Real FFT Using DSP-32C -----*/
init_dsp("ddj_32c.32c",&flag,&bs,np); /* download dsp code & init */
dl_float(input1,np,flag,b1); /* download float array */
set_dspbuf("o1", b1); /* set dsp buffer address */
dsp_dl_int(flag,IREALFT_C); /* execute "inv_realft" on dsp */
wait4dsp(flag); /* wait for dsp to finish */
ul_float(output1,np,flag,b1); /* upload results */
break;
case 8 : /*----- Calculate LOG(10)-[MAGNITUDE] using DSP-32C ----*/
init_dsp("ddj_32c.32c",&flag,&bs,np); /* download dsp code & init */
dl_float(input1,np,flag,b1); /* download float array */
set_dspbuf("i1", b1); /* set dsp buffer address */
set_dspbuf("o1", b2); /* set dsp buffer address */
dsp_dl_int(flag,LOGMAG_C); /* execute "inv_realft" on dsp */
wait4dsp(flag); /* wait for dsp to finish */
ul_float(output1,np,flag,b2); /* upload results */
break;
case 9 : /*------- Forward Real FFT Using DSP-32C App-Lib --------*/
init_dsp("ddj_32s.32c",&flag,&bs,np); /* download dsp code & init */
dl_float(input1,np,flag,b1); /* download float array */
set_dspbuf("o1", b1); /* set dsp buffer address */
dsp_dl_int(get_addr("stages"),ilog2(np));
/* download int */
dsp_dl_int(flag,RFFT_S); /* execute "rfft" on dsp */
wait4dsp(flag); /* wait for dsp to finish */
ul_float(output1,np,flag,b1); /* upload results */
break;
case 10: /*-------- Download Data To DSP-32C ---------------------*/
init_dsp("ddj_32s.32c",&flag,&bs,np); /* download dsp code & init */
dl_float(input1,np,flag,b1); /* download data from pc to dsp */
break;
case 11: /*--------- Upload Data From DSP-32C ---------------------*/
init_dsp("ddj_32s.32c",&flag,&bs,np); /* download dsp code & init */
ul_float(output1,np,flag,b1); /* upload results */
break;
default : set_err_return(BAD_FUNC_NUM); break;
}
}
/*--- This function returns the integer part of the log base 2 of x. ---*/
int ilog2(unsigned int x)
{
return( x >> 1 ? 1 + ilog2(x >> 1) : 0);
}
/*--- This function scales np elements of data1[] by scale_val. ---*/
void scale_data(float far * data1, float scale_val, long np)
{ long i;
for (i = 0; i < np; i++) {
data1[i] *= scale_val;
}
}
/*--- Function generates cosine data: data[i] = A cos((2 pi f/s) i) + d ---*/
void synth_cos(float far * data1, long np, float a, float d, float f, float s)
{ long i;
float theta, angle_step;
angle_step = 2.0 * PI * f / s;
theta = 0.0;
for (i = 0; i < np; i++) {
data1[i] = (a * cos(theta)) + d;
theta += angle_step;
}
}
/*--- log_mag ---*/
void log_mag(float far * i1, float far * o1, long bs)
{ long i;
long n;
n = bs >> 1;
o1[0] = log10(SQR(i1[0]));
for (i = 1; i < n; i++) {
o1[i] = log10(SQR(i1[2*i]) + SQR(i1[2*i+1]));
}
for (i = n; i < bs; i++) {
o1[i] = 0.0;
}
}
[LISTING TWO]
/*-----------------DDJ_MP87.C-------------------------------------------
| MathPak 87 FFT Function Group Execution Source File
| The makefile is: ddj_mp87.exe: ddj_mp87.c
| cl /AL ddj_mp87.c -link hqmcil.lib mpak87.lib
+----------------------------------------------------------------------*/
#include <hq_ci_.h> /* include function prototypes */
#include <mpak87.h> /* include math pak header */
#define BAD_BLOCK_SIZE -1 /* user defined error codes */
#define BAD_FUNC_NUM -2
int fft_stages(long fft_size); /*function prototype*/
void main(int argc, char *argv[])
{
int m; /* number of fft stages */
far_array_of_double o1; /* data pointer */
init_intfc(argc, argv); /* MUST do this before other interface functions */
o1 = get_data_out_ptr(1); /* get address of data */
if ((m = fft_stages(get_block_size())) == 0)
set_err_return(BAD_BLOCK_SIZE); /* won't happen if .fnc file correct */
else
switch(get_func_num()) {
case 1: rvfft(o1, m); break; /* real fft from MathPak library */
case 2: irvfft(o1, m); break; /* inverse real fft from MathPak */
default:
set_err_return(BAD_FUNC_NUM); /* won't happen if .fnc file correct */
break;
}
}
/* return the log(base 2) of the input, or 0 if input is not a power of 2 */
int fft_stages(long fft_size)
{ int rtn;
int sw_fft;
sw_fft = fft_size;
switch(sw_fft) {
case 8 : rtn = 3; break;
case 16 : rtn = 4; break;
case 32 : rtn = 5; break;
case 64 : rtn = 6; break;
case 128 : rtn = 7; break;
case 256 : rtn = 8; break;
case 512 : rtn = 9; break;
case 1024 : rtn = 10; break;
case 2048 : rtn = 11; break;
case 4096 : rtn = 12; break;
case 8192 : rtn = 13; break;
default : rtn = 0; break;
}
return(rtn);
}
[LISTING THREE]
/*--------------------DDJ_32C.C--------------------------------------
| This file is will run on the DSP add-in board after compilation
| by the AT&T C compiler d3cc. The makefile is:
| .c.o:
| d3cc -c $<
| .s.o:
| d3as -W -Q $<
| ddj_32c.32c: ddj_32c.o startup.o memory.map four1.o realft.o
| d3cc -lm -lap -o ddj_32c.32c ddj_32c.o four1.o realft.o \
| -s startup.o -m memory.map
+--------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#include <nr.h>
#define PI 3.1415926
#define SQR(a) ((a)*(a))
asm(".global i1, i2, o1, o2");
asm(".global funcnum, bs");
asm(".global amp, dco, freq, samprate");
short funcnum;
short bs;
float *i1, *i2, *o1, *o2;
float amp, dco, freq, samprate;
/*------------------------------------------------------------*/
short fft_stages(fft_size)
short fft_size;
{
short rtn;
switch(fft_size)
{
case 32 : rtn = 5; break;
case 64 : rtn = 6; break;
case 128 : rtn = 7; break;
case 256 : rtn = 8; break;
case 512 : rtn = 9; break;
case 1024 : rtn = 10; break;
case 2048 : rtn = 11; break;
case 4096 : rtn = 12; break;
case 8192 : rtn = 13; break;
default : rtn = 0; break;
}
return(rtn);
}
/*------------------------------------------------------------*/
main()
{ register float scal;
short n;
float *data1, *data2, temp;
register short i,j;
while (1) {
funcnum=0;
while (!(funcnum))
; /* Wait for PC to Download Function Number */
n = bs >> 1;
switch(funcnum) {
case 1: /*--------- Convert to DSP Format ------------*/
dsp32(bs,o1); /*IEEE-->DSP */
break;
case 2: /*--------- Convert to IEEE format -----------*/
ieee32(bs,o1); /*DSP-->IEEE*/
break;
case 3: /*--------- Synthesize Cosine ----------------*/
scal = 2.0 * PI * freq / samprate;
j = 0;
data1 = o1;
for (i=bs; i-- > 0; j++) {
*data1++ = (amp * cos(scal * j)) + dco;
}
break;
case 4: /*------- Forward FFT ------------------------*/
data1 = o1;
data1--;
realft(data1,n,1); /*from numerical recipes*/
break;
case 5: /*------- Inverse FFT ------------------------*/
data1 = o1;
data1--;
realft(data1,n,-1); /*from numerical recipes*/
/* scale by 1/n to retain original amplitude*/
data1 = o1;
scal = 1.0 / n;
for (i=bs; i-- > 0; data1++) {
*data1 = *data1 * scal;
}
break;
case 6: /*----- Calc LOG-MAGNITUDE (data output from NR-realft)--*/
o1[0] = log10(SQR(i1[0]));
temp = log10(SQR(i1[1]));
for (i=1;i<n;i++) {
o1[i]=log10(SQR(i1[2*i])+SQR(i1[2*i+1]));
}
o1[n] = temp;
for (i=n+1; i<bs; i++) {
o1[i] = 0.0;
}
break;
case 7: /*------ Forward FFT ----------------------------------*/
data1 = o1;
rffta(bs,fft_stages(bs),data1); /*uses AT&T app lib*/
break;
default : break;
}
}
}
[LISTING FOUR]
/*------ file DDJ_32S.S-----------------------------------------------
| Assembly language version of FFT test. The makefile is as follows:
| ddj_32s.32c: ddj_32s.s
| d3make -o ddj_32s.32c -M6 -Q -W ddj_32s.s
+---------------------------------------------------------------------*/
#include <dspregs.h>
/*--------------------------------------------------------------------*/
.global i1,o1,o2
.global funcnum, bs, stages
.global endofcode
/*--------------------------------------- initialization -------------*/
r5 = 0x0003
pcw = r5
ioc = 0x30CC0
/*--------------------------------------- wait until funcnum != 0 ----*/
begin:
r5 = *funcnum
nop
if (eq) goto begin
nop
/*--------------- switch on funcnum, which is set to function number --*/
/* func 1: IEEE-->DSP */
r5 = r5 - 1
if (eq) goto do_dsp32
/* func 2: DSP-->IEEE */
r5 = r5 - 1
if (eq) goto do_ieee32
/* func 3: invoke rffta, from AT&T app library */
r5 = r5 - 1
if (eq) goto rffta
/* func 4: calc magnitude squared */
r5 = r5 - 1
if (eq) goto do_mag
nop
/* func 5: calc log */
r5 = r5 - 1
if (eq) goto do_log10
nop
/* illegal function number */
goto finished
nop
/*---------------------------------call to _rffta -----------------*/
rffta:
r2e = *o1
r1 = *bs
r3 = *stages
*fft_b = r2e
*fft_n = r1e
*fft_m = r3e
call _rffta (r14)
nop
fft_lv: int24 localV
fft_n: int24 0
fft_m: int24 0
fft_b: int24 0
.align 4
goto finished
nop
/*---------------------------------calc magnitude------------------*/
do_mag:
r8 = *bs
r10e = *i1
r8 = r8 - 2
r11e = *o1
nop
a0 = *r10++ /* DC value */
nop
a2 = *r10++ /* Nyquist */
*r11++ = a0 = a0*a0 /* save DC mag */
magloop:
a0 = *r10++
nop
a1 = *r10++
a0 = a0*a0
nop
a1 = a1*a1
nop
nop
*r11++ = a1 = a0 + a1
if (r8-->=0) goto magloop
nop
*r11++ = a0 = a2*a2 /* save Nyquist mag */
goto finished
nop
/*---------------------------------calc log10------------------*/
do_log10:
r12e = *i1
r11e = *o1
r13 = *bs
r10e = in
r9e = out
r13 = r13 - 2
logloop:
*r10 = a3 = *r12++
call _log10 (r14)
nop
int24 localV
int24 in, out
.align 4
*r11++ = a3 = *r9
if (r13-->=0) goto logloop
nop
goto finished
nop
/*--------------------------------------call _ieee32 converter-----*/
do_ieee32:
r1e = *o1
r2 = *bs
*o1_ieee32 = r1e
*bs_ieee32 = r2e
call _ieee32 (r14)
nop
int24 localV
bs_ieee32: int24 0
o1_ieee32: int24 0
.align 4
goto finished
nop
/*--------------------------------------call _dsp32 converter-----*/
do_dsp32:
r1e = *o1
r2 = *bs
*o1_dsp32 = r1e
*bs_dsp32 = r2e
call _dsp32 (r14)
nop
bs_dsp32: int24 0
o1_dsp32: int24 0
.align 4
goto finished
nop
/*-------------------------------finished, set funcnum=0 -----------*/
finished:
r1=0
goto begin
*funcnum=r1
bs: int 256
stages: int 8
funcnum: int 0
.align 4
i1: int24 0x2000
o1: int24 0x2000
o2: int24 0x3000
.align 4
localV: 2*float 0.0
in: float 0.0
out: float 0.0
max: float 0.0
scalefac: float 0.0
#include <_rffta.asm>
#include <_log10.asm>
#include <_ieee32.asm>
#include <_dsp32.asm>
/*------------------------------------- mark end of code------------*/
endofcode: int 0xDEAD, 0xC0DE